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Abstract 

Wc consider a point process i + ^i, where j e Z and the ^^'s are i.i.d. random vari- 
ables with variance a^. This process, with a suitable rescaling of the distribution of 
i^i's, converges to the Poisson process in total variation for large a. We then study a 
simple queueing system with our process as arrival process, and we provide a complete 
analytical description of the system. Although the arrival process is very similar to 
the Poisson process, due to negative autocorrelation the resulting queue is very differ- 
ent from the Poisson case. We found interesting connections of this model with the 
statistical mechanics of Fermi particles. This model is motivated by air traffic systems. 

Keywords: Queueing system, air-traffic congestion, non Poissonian arrivals. 
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1 Introduction 



The main aim of this paper is to define a stochastic point process to model the arrivals to 
a queueing system, and to compare its features to the Poisson process. 
It is well known that the memoryless property of the Poisson process simplifies many tech- 
nical steps in the analysis of queueing systems, but there are arrival processes where such 
an assumption is not completely satisfied. In particular, we have in mind air traffic models. 
In recent times the dramatic increase of air traffic stimulated a large number of studies 
concerning the optimization of congestion management. From the point of view of classical 
queueing theory the system is difficult to study, mainly because it is hard even to define 
the basic quantities of the theory. For instance it is clear that there is some congestion for 
landing aircrafts, since they have to follow some holding paths, but it is not easy to quantify 
the actual time spent in queue or even its instant length. On the other hand, even assuming 
that the parameters of the system are known, it is not clear what kind of point processes 
are suitable to describe arrivals and service times. A common hypothesis in literature is 
to assume that arrivals are very well modeled by a Poisson process. This assumption, to 
our knowledge, goes back to the 70's when Dunlay and Horonjeff gave in [7j a number of 
theoretical and statistical arguments to justify the Poissonian hypothesis, and , since then, 
several other statistical studies have supported the same results. Even recently, see [6], a 
very careful study of the interarrival times of aircrafts to major US airports shows a small 
difference between the Poisson and the observed distribution, i.e. the actual arrivals are 
slightly less random than Poissonian ones, but the difference is quite small in all observed 
airports. On this ground, in various papers, see for instance [8], [9] and [lOj and reference 
therein, Poisson arrivals have been assumed in the analysis of judicious management of 
service times. It should be stressed that in all these papers the statistical validation of 
the Poissonian hypothesis has been based on computations on time scales smaller than the 
intrinsic randomness of the system. 

Stochastic models of aircraft arrivals based on statistical analysis and on simulations 
have a long history. As a first attempt, Barnett et al. [1] studied the arrivals to Boston 
Logan Airport. A version of the alternative model of arrivals we propose in this paper was 
introduced and studied numerically in [Ij. The model is refined in [3], where seasonal and 
daily effects are taken into account to describe random delays of departure times and, with 
these corrections, the model is quite accurate in its predictions. The key feature of the 
model is a soft a-priori scheduling of arrivals: indeed, both in US and in Europe, aircrafts 
are supposed to take off and to land by a schedule dictated by the capacity constraint of 
the runways, and by the assumption that each aircraft would land in a very narrow time 
slot. However, on the day of operations, an aircraft will be declared "on time" if it lands in 
a time interval larger than ten times the original slot. In this sense the scheduling should 
be considered "soft". The fact that arrivals are prescheduled clearly makes the Poissonian 
hypothesis questionable, but this is usually neglected, on the basis of the statistical studies 
mentioned above. However the predictions of the queueing theory give in general very rough 
estimates of the actual queue length. Moreover if we forecast a reduction of the intrinsic 
variability of arrival times, which could be achieved by various technical improvements (e.g. 
a rescheduling closer to the actual arrival times, or an en-route control of the paths of the 
aircrafts), we can not use Poissonian arrivals to describe the system, because they depend 
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only on a single parameter p. 

The process we study below is an arrivals model with two features. First, it shows a 
pattern of arrivals very close to a Poisson process when we look at time scales smaller than 
the standard deviation of aircraft delays, second, it provides the distribution of arrivals on 
time scales larger or comparable to the standard deviation of aircraft delays. 

Thus, the aim of this paper is an attempt to study more rigorously the features of arrival 
process presented in [4], which we suitably generalize, and to understand its analytical 
properties. 

Moreover, we show, both analytically and numerically, that the congestion related to 
this process is very different from the congestion of a Poisson process, on any time scale. 
This is due to the negative autocorrelation of the process, as we prove explicitly. It is worth 
to outline that the queueing models with Poisson arrivals have in general probabilities to 
have n customers in the queue that decay much slower than the probabilities observed in 
the air traffic. Our model gives a tail of the distribution much thinner, and more similar to 
the observed data. 

The analytical description of the system clarifies many interesting features of this kind 
of traffic: for heavy traffic the system has a long memory of the initial conditions; its 
description is obtained by the superposition of two processes, living on different time scales. 
This give the possibility to investigate also systems with slowly variable traffic intensities. 

The paper is organized as follows: in section 2 we describe our arrival process, and we list 
some results on the comparison to the Poisson process. In section 3 we present a simplified 
computation, obtained neglecting the autocorrelation of the process. The congestion levels 
according to this approximated process, assuming deterministic service (landing) times 
and a single server (runway), are quite different from the congestion according to Poisson 
arrivals. However we show numerically that such approximation is bad for very congested 
systems, where the actual level of congestion is not well described if the autocorrelation of 
the process is neglected. In section 4 we describe completely our queueing system at the 
price to enlarge suitably the state space of the Markov chain describing it. It turns out that 
for our process we have a finite value of the expected queue length even in the critical case 
^ = 1, while the Poisson queue diverges. Starting from the results on the critical case, we 
propose an approximation scheme that works very well for highly congested {q near to 1) 
systems. In this description a nice connection with the statistical mechanics of Fermi gas 
emerges quite naturally. Section 5 is devoted to conclusions and open problems. 

2 Description of the model: the arrival process 

In this section we want to introduce an arrival process, which we will call pre- scheduled 
random arrivals (PSRA) process, and to study its main features. The PSRA process is 
defined as follows. Let be the expected interarrival time between two clients, we define 
G M the actual arrival time of the z-th client by 



ti = T + ii 



i G Z 




where ^j's are i.i.d. random variables. 
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If the ^i's are uniform, the model is the actual arrival times process introduced in [3] 
without cancellations and pop-ups. We will show later that cancellations and pop-ups 
can be easily integrated into the process. From now on, we will assume that .^j's have 
continuous probability density f^'^\t) with variance cj^, and we will set without loss of 
generality = 0, since 7^ affects only the initial configuration of the system. 

The main aim of this section is to compare the features of the PSRA process to the Poisson 
process when a is large. It is well known, e.g. [21 p. 447], that the Poisson arrival process 
is defined by the fact that probabilities Pjj+i{At) = P{n{t + At) = j + = j) of a 

"jump" from the state j to the state j + 1 in the time interval (t, t + At] have the form 



Pjj+i{At) = P+{At) = XAt + o{At) 



(2.2) 



where A is a constant independent of t and j; A has the meaning of velocity of arrivals, i.e. 
denoting with ta the interarrival time, E{ta) = \- For PSRA the probability P{i,t,At) 
that the z-th client arrives in the time interval (t, t + At] is given by 



(2.3) 



&\x)dx (2.4) 



P(i, t,At) = P [t < - + <t + At 



P{t-^<ii<t + At-^ 



and, for small At, it may be written as 



P(i, t. At) = ^""^ ( t - ^ ) At + o(At) 



(2.5) 



By (j2.5p . the probability P~^{t, At) of a single PSRA arrival in the interval (t, t + At] is 
P+(t,At) = ^P(i,t,At)J];(l-P(i,t,At)) = 



E 



if {t-j-]At + c(At) 




/fMt-^) At + o(At) 



exp 



E 



f(^) it-±\At + o{At) 



A 



(2.6) 



Hence up to the first order in At the rate of arrival A(t) of the pre-scheduled random arrivals 
is defined by 



(2.7) 



This rate A(t) is periodic in t with period j. However we are interested in the dependence 
of A(t) on a, in particular when a is large with respect to j. To prove limit properties for 
our process, we have to specify the way we want to send a to infinity. We will require the 
following scaling property for the density f^'^\t). 
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Assumption 2.1. The probability density of ^ has the form 

/f)(t,a2) = i/^(t/a) (2i 



i.e. it is the rescaling of a well defined continuous density f^{t) with finite variance. We 
will also write maxfgjR f^{t) = M . 

This assumption is introduced in order to exclude pathological ways to send a to infinity, 
as, for instance, to consider a bimodal distribution with fixed maxima, see figure [TJ 





Figure 1: A bimodal distribution with fixed shapes shifting to infinity for o" — > oo. 
By our assumption, it follows that, in the limit cr very large the expression 

1 Jrr) ( . i 



i?(a,l/A):= (2.9) 



is the Riemann integral of the function f^'^\t). 
For example, let ^ be Gaussian N{0,a^), 



where Xi = and Ax = ^ and the limit is for cr —> oo. 

For any random variable rescaled in the above sense it is clear that the result 

lim i?(cr, 1/A) = 1 (2.10) 

cr— >oo 

holds, and therefore, in the same limit, 

lim A(t) = lim Ai?(o-, 1/A) = A (2.11) 

cr— >oo cr— ►oo 

It is interesting, for Gaussian ^, to check numerically how fast the limit is reached. Tabled] 
shows it. For simplicity, we set A = 1. 

The graph in figure [2] shows that, in terms of rate of arrivals, the pre-scheduled random 
arrivals approach the Poisson process when a is suitably large. In particular for Gaussian 
variables with standard deviation a of order 1/A or more we have that A(t) is constant up 
to 6 digits. Note that for applications mentioned in the introduction, we do expect the 
standard deviation to be much larger than 1/A. Note also that the explicit structure of the 
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a 


A(0) 


A(O.l) 


A(0.2) 


A(0.3) 


A(0.4) 


A(0.5) 


.2 


1.994726 


1.760407 


1.210523 


0.651951 


0.292114 


0.175283 


.3 


1.340089 


1.274318 


1.103259 


0.894087 


0.726696 


0.663191 


.4 


1.085005 


1.068767 


1.026261 


0.973729 


0.931237 


0.915008 


.5 


1.014384 


1.011637 


1.004445 


0.995555 


0.988363 


0.985616 


.6 


1.00164 


1.001327 


1.000507 


0.999493 


0.998673 


0.99836 


.7 


1.000126 


1.000102 


1.000039 


0.999961 


0.999898 


0.999874 


.8 


1.000007 


1.000005 


1.000002 


0.999998 


0.999995 


0.999993 


.9 


1. 


1. 


1. 


1. 


1. 


1. 


1. 


1. 


1. 


1. 


1. 


1. 


1. 



a 


A(0.6) 


A(0.7) 


A(0.8) 


A(0.9) 


A(l) 


.2 


0.292114 


0.651951 


1.210523 


1.760407 


1.994726 


.3 


0.726696 


0.894087 


1.103259 


1.274318 


1.340089 


.4 


0.931237 


0.973729 


1.026261 


1.068767 


1.085005 


.5 


0.988363 


0.995555 


1.004445 


1.011637 


1.014384 


.6 


0.998673 


0.999493 


1.000507 


1.001327 


1.00164 


.7 


0.999898 


0.999961 


1.000039 


1.000102 


1.000126 


.8 


0.999995 


0.999998 


1.000002 


1.000005 


1.000007 


.9 


1. 


1. 


1. 


1. 


1. 


1. 


1. 


1. 


1. 


1. 


1. 



Table 1: 



density of ^ does not play any particular role, and similar results may be obtained with 
different distributions. However it is clear that a small dependence on t is always present in 
the expression of A(t), and hence it is difficult to obtain a quantitative comparison between 
the pre-scheduled random arrivals and the Poisson process on this basis. Hence we look at 
the distribution of the random variable n{t, t + T), number of arrivals in the finite interval 
(t, t + T]. Let us call Pi{t, t + T) the probability that the i-th client arrives in the interval 
{t,t + T]. Clearly 

Pi{t, t + T) = {x - '-^ dx (2.12) 

Given the probabilities Pi{t,t + T) we can write the generating function of the random 
variable n{t, t + T), and, defining q^n^ = P{n{t, t + T) = n) we get 

1n^= E llP^it,t + T)llil-pjit,t + T)) (2.13) 

I={il,...,in} i<^I 30 

where the sum runs over all the possible distinct subsets / of indices of cardinality n. By 
mean of this expression one obtains the generating function 

g('^)(z) = q^:h- = n(l + {z- l)p^{t, t + T)) (2.14) 

n>0 ieZ 

To take into account also the possibility of random independent deletion as in [3], let us 
outline here that a similar generating function can be introduced also when each arrival has 
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Figure 2: Behavior of the function X{a, t) 



an independent probabihty 1 — 7 to be deleted, and the complementary probability 7 to be 
an actual arrival. In other words, we construct the PSRA process for i G Z and then for 
each i we cancel the corresponding i-th arrival with independent probability 1 — 7. It is 
obvious that in this case the generating function is 

q^^-\z) = 9^^" = 11(1 + - * + ^)) (2.15) 

n>0 iGZ 

The expressions ()2.14p . (j2.15p are exact, they give us all the information on the distribution 
of n{t,t + T), and they depend explicitly on t and T. However we want to study q^"\z) and 
q!^\z) for large a, in the sense of the rescaling defined above, showing that they converge 
to a Poisson distribution with parameter AT and 7AT respectively. The main idea is to 
exploit the fact that, for large a, Pi{t, t + T) goes to zero as ^. 
We now prove the following results. 



Lemma 2.2. 



maocpiit, t + T)< ^ (2.16) 

i a 



Proof. 

p,{t,t + T) = J^ /W(^^__jdx = y^ f'f\s)ds = -j^ h[-)ds (2.17) 

by the Intermediate Value Theorem 

p,(t,t + T) = -f^l^]T< (2.18) 
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where 



Si 1' I 



Now we will use lemma [2^ to bound the generating function 



q^"^ {z } = exp 



J^ln(l + (z-l)p,(t,t + r)) 



□ 



(2.19) 



exp 



{z-l)Y^p,{t,t + T)(l + {z-l)pi{t,t + T) [ ds- 
\ Jo < 



(i + (z-i)(i-s)Pi(t,f + r))2 

(2.20) 

Lemma 2.3. With pi{t,t + T) defined as above, the sum in (|2.20p converges to XT 



lim y Pi(t, t + T)(l + {z- l)pi{t, t + T) [ ds 

a^oo ^-^ \ If, 



{i + {z-i){i-s)pi{t,t + T)y 



Proof. First we prove that 



lim y^pi{t,t + T) = XT. 



= AT 
(2.21) 

(2.22) 



Let us define T := ^ where K G Z"*" and < AT < 1. Then we can write 



t+ 



K-i I AT 



&\s)ds 



i I K-i I AT 



&\s)ds 



(2.23) 



The first term on the right hand side of (j2.23p is K. Let i = mK + Z, where / € Z+ and 
m G Z, 



E 



(171-1)^+; 



K-l 



The second term on the right hand side of (j2.23p converges to AT for o" — > co: 



E 



t+ 



K-i 



iez ^+A 



A'^ A / S \ 

/^(-jds (2.25) 



and, by the Intermediate Value Theorem we get 

^ I K — i I AT , 



0- ./t+£^ 



a Xa 
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as o" ^ oo, where the sum on the last equahty is the Riemann sum of f^{t). This ends the 
proof of (j2.22p . In order to complete the lemma we need to show that, uniformly in i, 



lim {z-l)pi{t,t + T) [ ds — 

a^oo Jq (1 



+ {z-l){l-s)pi{t,t + T))2 
but this follows from lemma 12.21 and from the fact that 

s 



+ {z-l){l-s)piit,t + T)y 
for any Pi{t, t + T) < 1/2 and \z\ < 1. 



< C 



□ 



Lemma 2.4. Let q{z) = exp(XT(z — 1)) be the probability generating function of the Pois- 
son random variable ( with intensity XT, and q^yiz) = exp(7AT(z — 1)) be the probability 
generating function of the Poisson random variable C with intensity 'jXT, then 

lim = qiz); lim = qJz) (2.27) 

cr— >00 (T— >00 ' 

Proof. Follows immediately from lemma [2^ □ 

Theorem 2.5. If q^^^z) — > q{z), then Y,n= olln"^ ~ qn\ — > as a ^ oo. The same result 
holds for the arrivals with random deletions. 

Proof. The proof follows from the continuity theorem for probability generating function 
see Feller [H p.280]. □ 



Hence the PSRA process converges in distribution to the Poisson process in total variation 
norm, and the same is true for PSRA process with independent random deletions. 
In order to show that the process has negative autocorrelation, we will compute the expected 
value, the variance Var{n) of the number n of arrivals in a time slot (t,t + T], and the 
covariance Cov{ni,n2), where ni and n2 are the numbers of arrivals in {t,t + T] and {t + 
T,t + 2T], respectively. We present the explicit computations in the case of simple PSRA 
process, but the same results are true with obvious modifications for PSRA process with 
independent random deletions. 

Let Xi{ti £ (^5^ + ^]) be the characteristic function of the event "client i arrives in the 
interval {t,t + T]", so that E(xi) = Pi{t,t + T), then the expected number of arrivals in a 
time slot {t,t + T] is 

E(n)=Ef^Xi] =Y.E{x^) = Y,Pi{t,t + T) 

\ i / i i 
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and also 



= ^Mt,t + T) + Y, Mi, t + T)pj {t + T,t + 2T) 

i i^j 

= ^m{t,t + T)+ lY.p,{t,t + T)] -^{p,(t,t + T))^ 

i \ i / i 

Then the variance is: 

yor(n) =E(n2)-(E(n))2 = f+T)-^(pi(t, t+T))^ = Y^p^{t^t+T){l-pi{t,t+T)) 

i i i 

and we see again that Var{n) — > AT in the limit a — > oo. Finally, let us define Xi^^ '■— 
Xiiti e{t,t + T]) and xf ^ := XiiU e{t + T,t + 2T]) 



Y^p^{t,t + T)pj{t + T,t + 2T) 



= Y^Mt, t + T)pj{t + T,t + 2T) - Ep*(*' i + T)Mt + T,t + 2T) 

SO that 

Cot;(ni,n2) =E(nin2) -E(ni)E(n2) = -^Piit^t + T)pi{t + T,t + 2T) 

i 

A negative covariance means that ni and n2 are inversely correlated, as we should expect 
in our arrival model: a congested time slot should be followed or preceded by a slot with 
lower than expected arrivals. Moreover, this is a clear indication that the hypothesis of 
independence for ni and 722, numbers of arrivals in different time slots, is not correct, unless 
we are in the limit o" — > 00. 



3 Queueing systems with PSRA process: independence ap- 
proximation 

In this section we want to try to use the classical results of queueing theory for a system 
in which the arrivals are described in terms of our PSRA, there is a single server and 
the service time is deterministic. For the air traffic applications the deterministic service 
(landing) times are obviously an approximation, but neglecting the mix of aircrafts the 
actual landing times have a low variability. 
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a 


T 


g{a,0) 


g{a,OA) 


g{cj,0.2) 


^^(^,0.3)) 


g{a,OA)) 


.2 


.9 


0.808534 


0.808534 


0.850089 


0.907951 


0.954826 


.3 


.9 


0.868214 


0.868214 


0.88048 


0.900153 


0.919615 


.4 


.9 


0.892048 


0.892048 


0.895086 


0.900001 


0.904914 


.5 


.9 


0.898654 


0.898654 


0.899168 


0.9 


0.900832 


.6 


.9 


0.899847 


0.899847 


0.899905 


0.9 


0.900095 


.7 


.9 


0.899988 


0.899988 


0.899993 


0.9 


0.900007 


.8 


.9 


0.899999 


0.899999 


0.9 


0.9 


0.9 


.9 


.9 


0.9 


0.9 


0.9 


0.9 


0.9 


1. 


.9 


0.9 


0.9 


0.9 


0.9 


0.9 


a 


T 


g{a,0.5) 


^^(^,0.6) 


g{a,0.7) 


g{a,0.8) 


g{a,0.9) 


.2 


.9 


0.9786 


0.9786 


0.954826 


0.907951 


0.850089 


.3 


.9 


0.931537 


0.931537 


0.919615 


0.900153 


0.88048 


.4 


.9 


0.907951 


0.907951 


0.904914 


0.900001 


0.895086 


.5 


.9 


0.901346 


0.901346 


0.900832 


0.9 


0.899168 


.6 


.9 


0.900153 


0.900153 


0.900095 


0.9 


0.899905 


.7 


.9 


0.900012 


0.900012 


0.900007 


0.9 


0.899993 


.8 


.9 


0.900001 


0.900001 


0.9 


0.9 


0.9 


.9 


.9 


0.9 


0.9 


0.9 


0.9 


0.9 


1. 


.9 


0.9 


0.9 


0.9 


0.9 


0.9 



Table 2: 



In order to study our queueing process we set a service time T and we define the instant 
traffic intensity g{a,t) = E{n{t,t + T)). In fig. [3] and table [2] we report numerical results 
for the convergence of g{a, t) to AT, granted by lemma 12.31 For simplicity we consider ^ 
Gaussian, and A = 1. In this case g{cr,t) converges as soon as a gets close to 1. 

We want to compare the average queue size in M/D/1 queueing system (Poisson arrivals) 
with a queueing system in which the arrivals are described in terms of PSRA. To do this 
we have to recall some standard results in queueing theory. Assuming to have a probability 
Qn to have n arrivals in a service time slot, and assuming the variables n to be i.i.d, our 
system is described by the so-called discrete time GI/D/1 queueing model. 

It is well known, see e.g. [5], that the stationary probabilities for the discrete time 
GI/D/1 queueing model are given by 

Po = {Po + Pi)Qo 

(3.1) 

n+l 



Pn — PoQn + ^ PkQn-k+1 



k=l 

The corresponding generating function is 



Piz) = (3.2) 
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<3^0 XT 
>^^K p(<3,0.1) 

X p(0.2,t) 



Figure 3: Behavior of the function Q{a,t). On the x axis we have time t for Q{0.2,t) and 
standard deviation a for ^(fi, 0.1). 



In the case of Poisson arrivals with traffic intensity g, Q{z) = q{z) = exp(^>(z — 1)). Denoting 
by N the average queue size, after straightforward computations we get 



N 



g(2 - g) 
2(1 - e) 



(3.3) 



Consider now the PSRA process. In this case we can try to compute ()3.2p by means of the 
generating function ()2.14p . This is obviously an approximation, since for PSRA arrivals, 
as it has been shown in Section 2, the number of arrivals in subsequent time slots are not 
independent. 



However, neglecting the autocorrelation, we have that Q{z) 
N{a, t) the average queue size we find 



q^"\z), and denoting by 



Nia,t) 



2 E^&P^(t^ t + T)- {E^ezMt, t + r))^ - E^ezPKt, t + T) 



(3.4) 



For a large N{a,t) becomes independent of t, and it converges to N by (j2.22p . Table [3] 
shows that for Gaussian ^ and A = 1 the convergence is quite fast. 

The results obtained by the formulas above are an approximation, because we neglected 
the (negative) autocorrelations, and we have to see when this approximation is reliable. As 
a matter of fact the PSRA process is easy to implement for numerical simulation; hence 
we can compare the PSRA average queue size N{a, t) obtained by numerical simulations 
to ()3.4p and (j3.3p . In figure H] A^(cj, t) is plotted as a function of a, for different values of 
g = 0.5,0.7,0.9, and t = 0.5. The dotted straight lines represent N obtained by ()3.3p for 
different values of g. As we can see from the graph, values of N{a,t = 0.5) for fixed g given 
by (j3.4p are larger than the corresponding ones obtained by simulation. Moreover, this 
overestimate becomes very important when g increases. Hence, as it was easy to guess, the 
negative autocorrelation plays an important role in the system when the traffic intensity 
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a 


T 


iV(a,0) 


A^(cj,0.1) 


A^(cj,0.2) 


A^(cj,0.3) 


A^(cj,0.4) 


A^(cj,0.5) 


.1 


.9 


0.89105 


0.89105 


1.00493 


1.04024 


1.02267 


1.00905 


.2 


.9 


1.61425 


1.61425 


1.58187 


1.51872 


1.42902 


1.32201 


.3 


.9 


2.26812 


2.26812 


2.21399 


2.10656 


1.95949 


1.83453 


.4 


.9 


2.75253 


2.75253 


2.68673 


2.57205 


2.44587 


2.36133 


.5 


.9 


3.03548 


3.03548 


2.9955 


2.92993 


2.86327 


2.82151 


.6 


.9 


3.24502 


3.24502 


3.23019 


3.20614 


3.18205 


3.16714 


.7 


.9 


3.43207 


3.43207 


3.42809 


3.42165 


3.41521 


3.41123 


.8 


.9 


3.59488 


3.59488 


3.59405 


3.5927 


3.59134 


3.59051 


.9 


.9 


3.73131 


3.73131 


3.73117 


3.73094 


3.73071 


3.73056 


1. 


.9 


3.84462 


3.84462 


3.8446 


3.84457 


3.84454 


3.84452 



a 


T 


Af(cj,0.6) 


Af(CT,0.7) 


Af(cj,0.8) 


Af(cj,0.9) 


N{a, 1) 


.1 


.9 


1.00905 


1.02267 


1.04024 


1.00493 


0.89105 


.2 


.9 


1.32201 


1.42902 


1.51872 


1.58187 


1.61425 


.3 


.9 


1.83453 


1.95949 


2.10656 


2.21399 


2.26812 


.4 


.9 


2.36133 


2.44587 


2.57205 


2.68673 


2.75253 


.5 


.9 


2.82151 


2.86327 


2.92993 


2.9955 


3.03548 


.6 


.9 


3.16714 


3.18205 


3.20614 


3.23019 


3.24502 


.7 


.9 


3.41123 


3.41521 


3.42165 


3.42809 


3.43207 


.8 


.9 


3.59051 


3.59134 


3.5927 


3.59405 


3.59488 


.9 


.9 


3.73056 


3.73071 


3.73094 


3.73117 


3.73131 


1. 


.9 


3.84452 


3.84454 


3.84457 


3.8446 


3.84462 



Table 3: 



becomes large. For air traffic applications q near to the critical value ^> = 1 is the interesting 
case. 

4 Queueing systems with PSRA process: autocorrelated ar- 
rivals 

As it is clear from the results of the previous section, neglecting the autocorrelation the 
computed average queue length is grossly overestimated in the interesting cases. If we want 
to describe the system only by the length of the queue, the presence of autocorrelation 
implies the loss of Markov property. In this section we show that if we enlarge suitably the 
state space we may keep the Markov property, and describe completely the autocorrelation. 
With this description some interesting features of the system are clarified, but at the moment 
we are able to compute explicitly the quantities of interest with some approximations. Such 
approximations, however, turn out to give almost negligible errors. 

To simplify the analytical treatment of the system, we will consider from now on densities 
f'^\t) of the random i.i.d. variables that are compact support, i.e. such that f^\t) = 
for \t\ > L for some L < oo. We are setting A = 1, and we take L G N. This implies 
that at a certain discrete time j the i'th customer is certainly arrived to the system for 
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Figure 4: Behavior of the function N{a,0.5), for different values of g. Dotted lines refer 
to Poisson arrivals, continuous lines refer to approximation (j3.4p . dashed lines refer to 
simulations. The simulations are run for a time sufficiently long to have fluctuations on the 
result negligible in the scale of the figure. 

all i < j — L, while for all i > j + L it is certainly not yet arrived. Hence to completely 
describe the state of the system we have to specify, beside the number n of customers 
waiting in queue right before the service at time j is delivered, also a finite set Ij of i's, 
Ij C {j — L + j + L — 1}, that are the customers that are already arrived at the service at 
time j. Note that the customers in the set Ij are not necessarily already served at time j, or, 
in other words, the set Ij is the set of the customers with indices in {j — L + l,...,j + L — 1} 
that are in the queue at time j, or that are already served at time j. Note also that 
< \Ij\ <2L — 1. Finally, we want to outline that due to the independence of the ^'s /j+j 
is independent of Ij for all i > 2L. 

We will treat first the case g = 1, or in other words, the case A = T = 1 in (I2.14p . This 
special case is important for several reasons. First, we will prove that for PSRA arrivals 
the system has a finite average queue length, showing that, even if the PSRA process tends 
in distribution to the Poisson process, for finite variance of the ^'s the two systems are 
deeply different. Second, we will show that in the g = 1 case there is a conserved quantity 
in the system, when the stationary distribution is reached. Third, it is possible, using an 
interest interpretation of the system in terms of Fermi statistics, to compute the (very long) 
time needed to the system to reach the stationary distribution. Fourth, and maybe more 
important, on the basis of this computation it is possible to approximate efficiently the 
distribution of the length of the queue even for g < 1. 
Hence, we fix = 1 and we start from the obvious relation 

n(j + 1) = n(j) - (1 - 5„(,)o) + m(j) (4.1) 

where n{j) is the length of the queue immediately before the service at time j, m{j) is the 
number of customers arrived in the time slot [j,j + 1), and the term (1 — (5n(j)o) indicates 
the fact that if there is some customer in the queue at time j, i.e. n{j) > 0, the first of the 
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queue is served, while if n{j) = then n{j + 1) = m{j). 
Now we observe that with our notations we can write 



m(j) = |/,+i|-|/,-| + l (4.2) 

This relation can be shown as follows: the total number na[j) of customers arrived to 
the service from a certain fixed time, say from time 1, to time j, is obviously na{j) = 
j—L+\Ij\, because all the customers k up to customer j — L'th are already arrived, due to the 
compactness of the support of f^'^\t), while for k > j — L the number of arrived customers 
is \Ij\ by definition. Hence m{j) = na{j + 1) — na{j) = j + 1 — L + \Ij+i\ — j + L — \Ij\ = 
\Ij+i\ + Putting dlS]) intodiH) we obtain 

n{j + 1) = n(j) + - \Ij\ + 5„(j-)o (4.3) 

This relation shows that the quantity a{j) = n{j) — is constant during a busy period, 
and it increases by 1 at the end of each busy period. This implies that the stationary 
distribution is reached once a > 0. If the initial value of a is strictly positive, the value 
n[j) = is never realized, and then a remains constant and 

N = E{n) = a + E{\I\) (4.4) 

If the initial value of a is or it is negative, a sequence of busy periods is realized, giving 
in the end the value a = 1, and the expected queue length N = E{n) = 1 + £'(|/|). Once 
the stationary value of a > is reached, the probability distribution of n is given by 

= p{n = k) = P{\I\ = k-a) (4.5) 

giving the obvious result that k > a. The explicit expression of the Pk depends there- 
fore from the distribution of the |/|'s, and hence from the details of f^'^\t). This solves 
completely the stationary problem in the Q = 1 case. For application to the air traffic, 
however, it could be also interesting to study some non stationary features of the system: 
in particular we want to compute the probability to pass from some negative value of a to 
the following value q + 1. These quantities are interesting in this Q = 1 case because if the 
probability to reach the state n = for a given a < is much smaller that the inverse of 
the number of operation in a single day of traffic, it is very likely that the system remains 
on states n > 0. These probability to jump from a definite value of a to the following one 
are important also in the description of the ^ < 1 it will be explained below. 

Hence suppose that at time j the system is in the state n(j) = 0, with a given value of a < 0. 
Call t{a) the quantity such that n(j + i) > for all < i < t{a), and n(j + t{a)) = 0. 
t{a) is therefore the length of the busy period with starting value a. We are interested 
to the quantities T{a) = E{t{a)). By the definition of a we have that \Ij\ = —a + 1 
and that the instant j + t{a) is the first instant after j in which = —a, having 

> —a for all < i < t{a). To compute T{a) we should evaluate the probability 
P{\Ij^i\ = —a\\Ij\ = —a + 1). This probability are however hard to compute due to the 
conditioning. Here we introduce our approximation: we will measure T{a) in terms of 
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i.e. we neglect the conditioning. This approximation is reasonable for a such that P{\I\ = 
—a) <C in these cases we have to expect that the probability to have = 
— a||/j| = —a + 1) for i < 2L is very small, and since Ijj^i is independent of Ij for the 
greater values of i, that gives the bigger contribution to T(a), we have that the conditioning 
is almost ineffective. On the other side, for a such that P{\I\ = —a) > we have to expect 
a gross underestimate of = — a||/j| = —a + 1), and therefore a gross overestimate 

of T{a). We will return on this point later. 

We want now to compute explicitly P{\I\ = —a). We will write general formulas, valid 
for any density f^'^\t), and we will also consider a concrete probability distribution for the 

delays namely the case of f^^''\t) unifor m in [— L,L], in which many computations may 
be carried out explicitly. 

By straightforward computations one can see that 

^(|/|=0)= n (l-^cW) = 7^^e"'''^^ (4-7) 

i=-L+l ^ ' 



where the last approximation is valid for uniform ^'s, using Stirling formula, and 

1 

-L+l<ii<i2<...<ik<L- 



= P{\i\=^) E 7T^-7T^ (4.8) 

-L+l<ii<i2<...<ik<L-l 

where F^{t) is the probability distribution of the ^'s, and the last equality is again valid for 
uniform distribution. 

It is worthy to observe that (j4.8p may be interpreted as the canonical partition function 
of a Fermi system with 2L energy level and k particles, where the i-th level has energy 
log(Fg(i)) — log(l — F^{i)). With this respect many computational techniques may be used 
in order to compute the probabilities P{\I\ = k). Note that, in the approximation (j4.6p . 
once we are able to compute the quantities P{\I\ = k) we know also the expected values 
T(a). 

Let us list here a couple of possible way to evaluate P{\I\ = k) using the fact that, since 
it is possible to interpret it as a well known object in statistical mechanics, one can use 
computational results that are classical in that framework. The number of energy level, 
as mentioned above, is 2L. In real traffic context one should expect that this value is of 
the order 20 or 30. One of the available approximation of the quantity P{\I\ = k), i.e. 
the so called equivalence with the grand canonical ensemble, uses a method that is roughly 
speaking the Lagrange multipliers method, giving very good approximations for 2L large 
(see e.g. [HI chapter 5, section 53]). Since in our case 2L is not large enough to ensure the 
goodness of the approximation, it is much better to use an exact expression for P{\I\ = k), 
due to Ginibre. For completeness, and for the fact that it is quoted in a very implicit sense 
in [12] , we give the proof of this formula. 
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Calling Wi = j^^^, one can prove the following equality 

k I 



p{\i\ = k) = j2 E ^(^1' ■••'^■^) n 



/ = l<Ji<...<i; 



m=l i 



with 



cin,...,ji) = p{\i\=o) 



(-1) 



k-l 



ji jimi\...mk\ 



(4.9) 



(4.10) 



where rrii is the number of j's equal to i. To prove (j4.9p we observe that 



1 d!" 



P{\I\=k)=P{\I\=^)-—,J{{l+tw. 

i 

The quantity Y\^{1 + twj) can be expanded in series as follows 



t=o 



1 TT. 



t=0 



1 



E,E,=i(-i)^ 



-1 i-t^i)^ 

3 



t=0 



E 

i=0 



(-1) 



k-l 



E EE 

31. •••j'; m=l i 



t=0 
3m 



1 .^ (E-=i(-ir^fE. 



ID, 



«=1 



t=0 



E 

«=0 



(-1)'^ 



fc-Z 



E HE 

i<ji<...<j; m.=l j 



jrn mi\...mk\ 



which is (jOj) . 

We conclude then the discussion of the = 1 case observing that in a concrete framework 
of air traffic, if we want to avoid to have lost slot but we want to keep the queue as short 
as possible we have to choose initial condition in such a way that a is the smaller possible 
value such that T[a) > D, where D is the number of operations in a day. This value of a 
gives the corresponding value of the length of the queue using (|4.4p . 

A simple observation allows us to give an estimate of the average length of the queue also 
when Q < 1. Let us suppose that we impose the condition q < 1 keeping the time between 
two expected arrivals equal to the service time, but assuming that the arrivals are described 
by PSRA process with random deletion (see (|2.15p ). with probability of deletion equal to 
1 — ^. It is easy to realize that this corresponds to say that the value of a has a probability 
1 — Q to decrease by one. Hence we have this picture of our queueing system: the queue 
is described by a superposition of a slow varying process, the process that describes the 
value of a, and a fast varying process, the one describing the n for fixed a. If we are able 
to compute the distribution probabilities of the values of q, we can evaluate the expected 
length of the queue (and even its distribution) by (|4.4p . weighted with the probabilities of 
the various values of a. 
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In the unconditioned approximation (j4.6p . the computation of the stationary probabihties 
TTa of a is a standard task of the theory of the birth-and-death processes: the evolution of 
a is a discrete time birth-and-death process, with transition probabihties 

^ 1 — g = fia ifa' = a — 1 



P{\I\ = -a) = A„ if a' = a + 1 
1 — Xa — fia if a' = a 

otherwise 



and boundary conditions ^^^j^i = Aq = 0. We get the fohowing hnear system 

TT-L+l = 7r-L+l(l - A_L+l) + 7r_L+2^-L+2 

vTj = 7ri_iAi_i + vTj+i^i+i + 7rj(l - Xi - ^t) - L + 1 < i < 
TTo = 7r_iA_i + 7ro(l - no) 



whose solution is 

i 

k=-L+2 

The stationary distribution vr is defined by the normalization condition X^jVTj = 1, then 



This approximation is good for 1 — Q sufficiently small, because the probability to increase 
a = —L + 1 is much bigger than the probability to decrease it, and at the same time 
the unconditioned transition probabilities to increase a when a > —L + 1 are a good 
approximation of the actual transition probabilities. 

In the following figure we show the value of the expected length of the queue obtained by 
the formula 

N = J2^a{a + E^{\I\)) (4.12) 

a 

Note that Ea{\I\) is a-dependent, because in its computation we neglect the terms with 
|/| < —a, since they do not contribute to the evolution of the process with that value 
of a. As it can be seen from the figure, the estimate of the average length of the queue 
is extremely near to the simulations, also for highly congested systems. In the figure we 
have shown for completeness also the (wrong, for high g) values of the length of the queue 
computed by means of formula ()3.4p . which neglects the autocorrelations. 
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Figure 5: The length of the queue for highly congested systems, computed by means of 
numerical simulations (red line) and our analytical approximation (blue line). It can be 
seen that the uncorrelated approximation (black line) obtained by formula (j3.4p gives for 
these values q a gross overestimate. The simulations are run for a time sufficiently long 
to have fluctuations on the result negligible in the scale of the figure. 

5 Conclusions and open problems 

The main aim of this work is to study a stochastic process close to the Poisson process, but 
more suitable to describe the arrivals to a queueing systems when such arrivals are scheduled 
in advance, and some randomness is added to the schedule. We looked into this problem as 
an attempt to describe the congestion in air traffic systems, but the same construction can 
be used in different contexts. 

We found analytical results, in particular we showed that our process can be indistin- 
guishable from a Poisson process if one wants to study the distribution either of the number 
of arrivals or of the interarrival times in a time slot shorter than the standard deviation of 
the randomness imposed to the scheduled arrivals. 

However we have shown that from the point of view of the resulting congestion, due to the 
autocorrelation of this stochastic process, the queueing properties of this model are quite 
different from the analogous problem with Poisson arrivals. Interesting connection with 
the statistical mechanics emerged in the analytical solution of the problem. We proposed 
some approximation in our computations, but the results we obtained are in very good 
agreement with numerical simulations. An important question is the discussion of the 
accuracy of this description with respect to actual air traffic data. We have with this 
respect some preliminary results showing that the description of the distribution of the 
length of the queue using the PSRA as arrival process is much more accurate than the 
description assuming Poisson process, that is well known to be unfit. We hope that this 
study, that has to be fully developed in its computational aspects, may shed some light in 
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various unclear aspects of the air traffic modeling. 
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